*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*	This file generates Figure 2 using the MSE simulation output.
*	----------------------------------------------------------------------------
	local outname 		"${city}_RMSE"
*	----------------------------------------------------------------------------
	args ptype
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

	cap log close
	log using "${log}`sheetname'.smcl", replace

*	load

	clear
	set obs 1

		import delimited "${tables}${city}_sims_`ptype'_math1.csv", clear

		g vam = ""
		replace vam = "true" if _n == 1
		replace vam = "unc" if _n == 3
		replace vam = "conv" if _n == 5
		replace vam = "risk" if _n == 7
		replace vam = "rc" if _n == 9
		replace vam = "conv_3yrlag" if _n == 11
		replace vam = "unc_ivvam" if _n == 13
		replace vam = "conv_ivvam" if _n == 15
		replace vam = "risk_ivvam" if _n == 17
		replace vam = "rc_ivvam" if _n == 19
		replace vam = "conv_3yrlag_ivvam" if _n == 21

		keep if inlist(_n,1,3,5,7,9,11,13,15,17,19,21)

		ren v12 closure
		ren v25 mse
		ren v26 mse_var
		ren v27 mse_bias
		drop v? v1* v2* v3*

		g rmse = sqrt(mse)

		g bias_share = mse_bias/mse
		g bias = bias_share * rmse
		g var = rmse - (bias_share * rmse)

*	 plot

	g order = .
	replace order = 1 if strpos(vam,"unc")>0
	replace order = 2 if strpos(vam,"conv")>0
	replace order = 3 if strpos(vam,"risk")>0
	replace order = 4 if strpos(vam,"rc")>0
	replace order = 5 if strpos(vam,"conv_3yrlag")>0
	replace order = 6 if strpos(vam,"true")>0

	g barorder = .
	replace barorder = 4 if vam == "unc"
	replace barorder = 5 if vam == "unc_ivvam"
	replace barorder = 7 if vam == "conv"
	replace barorder = 8 if vam == "conv_ivvam"
	replace barorder = 10 if vam == "risk"
	replace barorder = 12 if vam == "risk_ivvam"
	replace barorder = 13 if vam == "rc"
	replace barorder = 14 if vam == "rc_ivvam"
	replace barorder = 15 if vam == "conv_3yrlag"
	replace barorder = 16 if vam == "conv_3yrlag_ivvam"

	gen ivvam = strpos(vam,"ivvam")>0

*	plot and save

	graph bar bias var if vam != "true", over(ivvam, gap(*.1) relabel(1  "OLS"  2  `" "IV" "VAM" "' 3 "{bf:OLS = IV}" ) label(labsize(small))) ///
		over(order, gap(*.3) relabel(1  "Uncontrolled" 2  `" "Conventional" "controls" "'  3 `" "Risk only" "controls" "' 4  `" "RC VAM" "controls" "' 5 `" "Conventional with" "older lagged scores" "')  sort(barorder) label(labsize(small))) stack nofill ///
		ylabel(,nogrid) legend(lab(2 "Variance") lab(1 "Bias") order(1 2) pos(6) r(1)) yti(Root mean squared error) ///
		bar(1,color(gs4) lwidth(none)) bar(2, color(gs9) lwidth(none))

	graph export "${figures}`outname'.pdf", replace
	graph export "${figures}`outname'.eps", replace

	log close
